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ABSTRACT 

We present a new calculation for the evolution of the 1-point Probability Distribution 
Function (PDF) of the cosmological density field based on an exact statistical treat- 
ment. Using the Chapman-Kolmogorov equation and second-order Eulerian pertur- 
bation theory we propagate the initial density distribution into the nonlinear regime. 
Our calculations yield the moment generating function, allowing a straightforward 
derivation of the skewness of the PDF to second order. We find a new dependency 
on the initial perturbation spectrum. We compare our results with other approxima- 
tions to the 1-pt PDF, and with N-body simulations. We find that our distribution 
accurately models the evolution of the 1-pt PDF of dark matter. 
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1 INTRODUCTION 

The standard scenario for the formation of structure in the 
Universe is via the gravitational amplification of primordial 
random Gaussian fluctuations generated in the Early Uni- 
verse during an Inflationary phase. An attractive feature of 
this scenario is its predictive power in determining the his- 
tory of mass perturbations from initial times up until the 
present day. In principle this should allow one to compare 
observation with theory. But while detailed predictions at 
high redshift for the Cosmic Microwave Background (CMB) 
are possible due to their linearity and dependence on well 
tested laboratory physics, predictions at lower redshift are 
complicated by nonlinear gravitational evolution and the 
physics of galaxy formation. One of the goals of cosmol- 
ogy is accounting for the statistical evolution of mass and 
galaxies up to the present day. 

In this paper we focus on the properties of the one- 
point density distribution function of matter. Interest in the 
distribution functions of cosmological fields has grown as 
they encode a great deal of information on initial condi- 
tions, gravitational evolution and galaxy bias. However the 
challenge is to separate out each of these effects. 

The evolution of the one-point distribution function 
has been computed by a number of authors using a vari- 
ety of techniques here we present a new approach, based 
on propagating the distribution function by the Chapman- 
Kolmogorov equation, incorporating nonlinear evolution to 
second order. In a further paper we shall develop the for- 
malism to include the effects of galaxy biasing and redshift- 
space distortions (Watts & Taylor 2000). 

The 1-pt PDF is a useful quantity in cosmology. 
While reasonably straightforward to measure (see Hamil- 



ton 1985, Alimi et al 1990, Szapudi, Szalay & Boschan 1992, 
Gaztafiaga 1992, 1994, Bouchet et al 1993, Szapudi, Meiksin 
& Nichol 1996, Kim & Strauss 1998) the 1-pt PDF embodies 
the full hierarchy of correlation functions and its measure- 
ment ensures that physical constraints such as the Lyapunov 
inequalities for the moments {{x"'^'') > {x''){x^), for any 
random variable x) are automatically satisfied. This is not 
necessarily the case for direct measurement of the moments. 
In addition the PDF offers a convenient way of dealing with 
Poisson sampling in the galaxy distribution. 

From a theoretical perspective, the 1-pt PDF is a conve- 
nient quantity to calculate if it is initially Gaussian. A num- 
ber of different methods have been developed to calculate its 
evolution. Fry (1985) first suggested calculating the proba- 
bility function by applying a hierarchical series suggested by 
the BBGKY equations to solve for the moment generating 
function. Bernardeau (1992) derived an exact expression for 
the evolution of the moment generating function. The gen- 
eralisation to top-hat filtered density and velocity fields is 
found in Bernardeau (1994), while Bernardeau (1996) stud- 
ied the 2-point cumulants and the PDF. However so far a 
generalisation of these results to include bias and redshift- 
space distortions in the PDF has not been presented. 

An approximation for the PDF can also be constructed 
from the first few cumulants, in the limit of small variance, 
by the Edgeworth expansion (Juszkiewicz et al., 1995 and 
Bernardeau & Kofman, 1995), where the cumulants have 
been derived directly from Eulerian perturbation theory 
(Bouchet et al., 1992) or Lagrangian perturbation theory 
(Bouchet et al. 1995). This can be extended to redshift- 
space using the Lagrangian perturbation calculations of 
Hivon et al. (1995) for the skewness. Colombi et al. (1997) 
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introduced an extended perturbation theory, where the re- 
sults of perturbation theory are allowed extra freedom and 
extrapolated to the nonlinear regime. 

As well as Eulerian perturbation theory other approxi- 
mations have also been applied. Kofman et al (1994) and 
Bernardeau & Kofman (1995) derived the PDF in the 
Zel'dovich approximation, while Hui, Kofman & Shandarin 
(1999) extended this to include the effects of redshift -space 
distortions. 

A more phenomenological approach was taken by Coles 
& Jones (1991) who approximated the properties of the PDF 
by a lognormal distribution while Colombi (1994) suggested 
an Edgeworth expansion about the lognormal distribution. 

In this letter we apply second-order perturbation the- 
ory to an initially Gaussian distributed density field to cal- 
culate the exact second-order characteristic function. We 
then numerically inverse Fourier transform this to yield the 
1-pt pdf. Our method therefore treats the propagation of 
probabilities exactly. 

The letter is organised as follows. In Section 2 we cal- 
culate the evolution of the cosmological probability distri- 
bution function and discuss the transformation to a discrete 
count distribution. In Section 3 we compare our results with 
the results of N-body simulations and with other distribu- 
tions in Section 4. Conclusions are presented in Section 5. 
We begin by describing second-order perturbation theory 
and calculating the probability distribution function. 



is the linear peculiar gravity field^ where V ^ is the inverse 
Laplacian. The trace-free tidal tensor is given by 



E,j{x,t) = V.VjV ^(5i(x,t) - -5x{x,t)5,j. 



(5) 



Equation (^) is the general result for an arbitrary cosmology 

wherefl k = D2/DI ^ -S/TJ^"^/*^^ (Bouchet et al 1992, see 
also Catelan et al. 1995), and 5\{x,t) — Di{t)ei{x), where 
Di{t) ^ (1 + z)~^ is the linear growth function (Peebles 
1980). 

2.2 The distribution function of initial fields 

In order to calculate the second-order density distribution 
function, P{S), we must first calculate the joint probability 
of each of the fields in equation P{Si, ri,g, E). Defining 
the parameter vector y — (5i, J?, g, E) the joint distribution 
function is given by the multivariate Gaussian 

^(y) = ((2vr)'^|detCi)^/^ H^^'y) ' ("^^ 
where 

C = ivy') (7) 

is the 12 X 12 covariance matrix that gives the degree of cor- 
relation between each of the components of y. The elements 
of the covariance matrix are 



2 THE COSMOLOGICAL PROBABILITY 
DISTRIBUTION FUNCTION 

2.1 Second—order perturbation theory 

Before shell-crossing, and in a spatially flat universe, the 
Eulerian density field, S{x, t), can be expanded in a series of 
separable functions; 



5{x,t) ^^5n{x,t) ^^D„{t)en{x), 



(1) 



where e„ is an n*''-order time-independent density field and 
Dn is an n"'-order universal growth function. It is a special 
property of this expansion that the spatial and time com- 
ponents are separable. In a universe with spatial curvature 
this expansion is only possible to second-order. 

The growth of perturbations in the nonlinear regime 
can be calculated by solving the continuity, Euler and Pois- 
son equations sequentially for each order in the perturbation 
expansion. To second-order the density field can be derived 
from linear quantities by the relation (Peebles 1980, Bouchet 
et al 1992) 



5 = 5i + -{2~K)5t 
where 



»7.g + i(l + K)ij2 



(2) 



■n{x,t) = VSi{x,t) (3) 
is the gradient of the linear density field and 



g{x,t) = -VV ^Si{x,t) 



(4) 



(EijEki) = i^'^o ('^'fc^J' + SiiSjk — ^SijSki^ , (8) 

with the remaining entries in the covariance matrix zero. 
These variances are defined as 



(9) 



where P{k) is the initial power spectrum of density pertur- 
bations. The components of this distribution are 



and 
detC 



^0 (1-7S) cr_i cri(T_i 



20 12 6 6 M 2n3 



3105»7r4 

The correlation parameter, 7^, is defined as 



(Ti cr_i 



(11) 



(12) 



providing a measure of the correlation between the linear ve- 
locity and density gradient fields. If we assume a power-law 
power spectra with a Gaussian cut-off, P{k) oc fc"e~* ^ , 



In this paper we define AirGpo = 3/2nH^ = 1 and the expan- 
sion parameter a{t) = 1, since our final distribution function will 
be dimensionless. 

t We define k following Bouchet et al. (1995) and Hivon et al. 
(1995), but note that this difi'ers from the definition of Bouchet 
et al. (1992) who define k ^ 2/14^-2/63 
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where n is the spectral index and R some arbitrary length 
scale, then 



7^ = 



n + 1 



(13) 



n + 3' 

with the constraint n > —1 to insure convergence of the 
velocity field. The correlation parameter must be positive 
definite since under gravity matter will be displaced from 
low- to high-density regions. The special case of n ^ — 1 
for a power-law initial spectra results in 71, = because the 
velocity field diverges on small scales. In fact the diverging 
velocities form an incoherent Gaussian random field which 
is uncorrelated with the density field (Taylor & Hamilton 
1996). 



2.3 Propagation of the density distribution 
function 




(1+6) 



The distribution function can be propagated to later times 
by the Chapman-Kolmogorov equation, 



P{x)= / dyWix\y)P{y) 



(14) 



where W{x\y) is the transition probability from a; to y. In 
the case of a deterministic process, such as the one consid- 
ered here, the transition probability reduces to a delta func- 
tion restricting the number of possible paths of evolution to 
one. This transition probability is given by 

W{5\y) =5o[5-5^- i(2 - k)5? + ^.g - 1(1 + k)E''] (15) 

Inserting this into equation ( p^ we find 

P{5)^{SD{5-S{y)))y. (16) 

where 5{y) is the right hand side of equation Hence for 
deterministic transitions the Chapman-Kolmogorov equa- 
tion becomes the expectation value of the delta function. 
This is a well-known result from probability theory (eg van 
Kampen 1992). 

Fourier transforming the delta function we find 



where 



dj Q{J) exp(iJ5) 



dSP{5) exp {-iJS) 



(17) 



= {exp(-iJ(S(i/)))j 



(18) 



is the characteristic function. In the second line we have 
used equation (^^ to write the characteristic function as 
an expectation over all the linear fields, y. This expression 
reduces to a set of multivariate Gaussian-type integrals that 
we can easily evaluate yielding 



g{J) = e(J)exp 



where 



J o-Q 



2(1 + iaiJ) 



(19) 



e(J) = {l+iaiJ)~''/^{l+ia2J)-''^^il-ia-iJ+aiJ^)-^^\20) 
and where the a coefficients are given by 



Qi = -(2 - K)ao, 



02 



15 



(1 + K)cro, 




Figure 1. The evolution of the 1-pt distribution function, P(<5) 
calculated from second-order perturbation theory. The variances 
are erg = 0.15, 0.20, 0.25, 0.3. In the lower plot we use logarithmic 
axis to emphasis the tails of the distribution. 



03 = -ctq, 



ill 

9 



yi) 4 

; — o-Q 



(21) 



The probability distribution can then be found by numer- 
ically integrating equation ([l^). Equations (|ig|), ( |2c| ) and 
( pl| ) are the central results of this letter. 

The characteristic function for S can be separated into 
characteristic functions for each term in the summation in 
equation (^). This is to be expected, since the distribution of 
a summation of random variables is equal to the convolution 
of their individual distributions, and so the characteristic 
functions just multiply. We expect this to be true for all or- 
ders of the perturbation series. But rather than converging 
to a Gaussian distribution, as suggested by a naive applica- 
tion of the Central Limit Theorem, the final distribution is 
driven away from it by the correlations induced by gravity. 

Figure 1 shows the evolved density distribution function 
for a range of ao . We use linear axis for the top plot empha- 
sising the peak and logarithmic axis in the lower plot empha- 
sising the tail of the distribution. As expected, the shape of 
the distribution is very nearly Gaussian when the variance 
is small, becoming very rapidly non-Gaussian for higher val- 
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ues. For high variances the probabihty density does not drop 
to zero at 5 = — 1, since in second-order perturbation theory 
the density field can be negative, generating non-vanishing 
regions with 5 < —1. This is also true of linear theory where 
there are always negative density regions for Gaussian initial 
conditions. This is a feature of Eulerian distribution func- 
tions. Those calculated in Lagrangian perturbation theory, 
or the lognormal, have positive definite densities at all times. 

In Figure ^ we demonstrate the effects of the correla- 
tion parameter, 7^, on the PDF. The main effects are an in- 
crease in volume of underdense regions with a corresponding 
decrease in extremely underdense regions when 7,^ is high. 
For low 7^, the underdensities are smaller and deeper. The 
physical reason for this is that the rj.g term in second-order 
perturbation theory deals with the evacuation of the voids, 
rather than the amplification of peaks. When 7„ is low this 
term is weakened and voids tend to be smaller and deeper, 
as they would be if the linear field were extrapolated. In- 
creasing this correlation widens the voids, but makes them 
shallower to help satisfy the 5 > —1 constraint. This effect 
is small for CDM-type initial power spectra where the cor- 
relation parameter is in the range 0.55 < < 0.65 over a 
wide range of scales. 



2.4 Skewness from the characteristic function 

Since our characteristic function is correct to second-order 
a useful check is to calculate the variance and skewness to 
compare with previous estimates. Taking the derivatives of 
the characteristic function with respect to J and setting J 
equal to zero yields the moments of the evolved density dis- 
tribution function, 

Qn 

(22) 



d[iJ]" 



g{j = o) = {5"). 



In the literature it is common to express the moments in the 
form of the moment parameters 



(23) 



where {5")c is the connected or irreducible part of {5"). The 
irreducible moments, or cumulants, can be generated by 



Qr, 



d[iJY 



■lne(J = 0) = {5"). 



(24) 



From this we can calculate the second and third order con- 
nected moments of our distribution function; 



(25) 



(26) 



52 = 1 
and 

53 = 2(2-.) 

to lowest order, reproducing the results of Peebles (1980). 
Hence our PDF leads to the correct second-order skewness. 
The intrinsic effects of the shape of the power spectrum via 
the correlation coefficient 7,, are to higher order. 



2.5 The discrete distribution function 

In reality the density field of galaxies is not a continuous 
function, but a discrete distribution. To account for this 
it is usual to assume Poisson sampling of the contirmous 





Figure 2. The variation in the peak (upper) and tail (lower) of 
the 1-point density distribution function as a function of the cor- 
relation parameter, 7^/. The values used were 7^ = 0.65 (solid 
line), 0.55 (dashed), 0.45 (dot-dashed) and 0.35 (dotted). The 
effect is small for CDM-type initial power spectra where the cor- 
relation parameter is in the range 0.55 < 7^ < 0.65 over a wide 
range of scales. 



density field as a crude approximation to galaxy formation 
processes. Any variation is attributed to biased galaxy for- 
mation, which modulates the underlying function. We shall 
consider this elsewhere (Watts & Taylor, 2000). 

The continuous density distribution function can be 
transformed to a discrete form by the expectation 



Pin) 



(27) 



where n is the mean galaxy count and the expectation is 
taken with respect to the nonlinear density distribution. Ex- 
panding this in terms of the characteristic function and tak- 
ing the expectation we find 



Pin) 



dJ 



GiJ) (1 



iJ 



exp {—iJ) 



(28) 



In the limit n 00 this returns the continuum distribution 
where S = n/n — 1. We note that again that the discrete 
distribution has a non- vanishing value at n = 0, P{n = 0), 
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since the probability of finding a no galaxies within a cell is 
finite. This is the Void Probability Function (White 1979). 

It is also useful to have the generating function for the 
discrete moments of this distribution, Gn(fc). Following some 
straightforward calculation we find that 



Gn{k) = G[in{e — 1)] exp n(e 



(29) 



The discrete moments can then be found by differentiating 
G„(fc); 



in 



d"'Gn{k = 0) 
d{ik)"^ 



(30) 



Again the connected moments can be found by differentiat- 
ing In G{k) with respect to ik. 



3 COMPARISON OF RESULTS WITH 
N-BODY SIMULATIONS 

In this section we show a comparison between our theo- 
retical PDF and the counts in cells PDF found from cos- 
mological n-body simulations. We used a version of Hugh 
Couchman's (1991) Adaptive P^M code altered by Peacock 
& Dodds (1994) to allow simulations of low density open 
and flat universes. 

The simulation volume was a periodic cube of comov- 
ing side 200/i~^Mpc containing 100'^ particles. We chose a 
CDM-type linear power spectrum of Gaussian initial per- 
turbations (Bardeen et al. 1986), normalised to match the 
observed abundance of clusters with linear variance given by 
as = 0.6^"° '^^ (Viana & Liddle 1996). The simulation was 
carried out on a 128^ Fourier mesh. 

To measure the PDF from the numerical simulations 
we smoothed the discrete particle distribution with a Gaus- 
sian filter of radius R. The PDF was then found from the 
smoothed density field evaluated on a 128^ grid. Since the 
effective radius of the binning grid was much smaller than 
that of the filter radius we expected negligible contribution 
to the smoothing from it. 

The choice of variance to use in our model is slightly 
ambiguous, given that we do not treat smoothing exactly. 
Hence it is better to match variances rather than smoothing 
scales. However the variance in second-order perturbation 
theory is the same as in linear theory so we could use either 
the linear input variance from the simulations, or take the 
measured nonlinear variance when making a comparison. In 
practice we found that all of the models provided a better 
fit if the nonlinear variances where used. This is a well know 
effect and is the basis of the extended perturbation theory 
approach of Colombi et al (1997). 

Figures (^, top) and bottom) show the results of 
numerical simulations (points) alongside our second-order 
PDF (solid line), the lognormal model (Coles & Jones 1991, 
dotted) and the Edgeworth expansion (Juszkiewicz et al. 
1995, Bernardeau & Kofman 1995; dashed). Each plot shows 
the PDF for two different variances, ctq — 0.2, smoothed 
on a scale of i? = 20/i~^Mpc (open circles) and 0.3 (filled 
squares) smoothed on 13.5/i~^Mpc. 

The effects of shot noise were taken into account as 
shown in Section 2.3 with the mean particle density n chosen 
to coincide with the mean particle count in a boxes with 
side I = \^R. Since n was a very large number in all cases. 





(1+6) 

Figure 3. Comparison of the second-order PDF (solid line) with 
the results of N-body simulations (points). The variances are 
(TO = 0.2 (open triangles) and 0.3 (filled squares). Also plotted 
are the corresponding distribution functions for the lognormal 
model (dotted) and Edgeworth expansion (dashed). 



shot noise effects were extremely small. This would not be 
the case in a real galaxy survey with lower number counts. 
The remaining factor controlling the shape of the theoretical 
PDF was 7^. For linear CDM power spectra the calculated 
7^ was approximately 0.55 for a wide range of smoothing 
radii and it was this value that was used to prepare all the 
plots in this section. 

Error bars on the N-body data points are the standard 
deviation over 5 independent simulations with identical cos- 
mological parameters. 

The linear axes of figures (|^, top) shows the accuracy 
around the peak of the distributions while the logarithmic 
scale of figure (^, bottom) shows the tails of the distribu- 
tions. At low variance (ctq = 0.2) we find very good agree- 
ment between all of the models and the simulations. Overall 
we found that the lognormal model fitted the simulations ex- 
tremely well for all values of the variance and around both 
the peak and the tails. As has been remarked elsewhere we 
regard this as something of a fiuke (Bernardeau & Kofman 
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based on the agreement we find between it and our N-body 
simulations. 

Over the range 0.4 < 5 < 1.5 all of the models agree to 
within a few percent, with the exception of the second-order 
Edgeworth expansion. But for densities below this regime, 
our second-order PDF overpredicts because the positive 
density condition is only weakly met. Both the Zel'dovich 
and Edgeworth models underpredict the probability of low 
density regions. Above the regime of good agreement our 
model again overpredicts, again because we over-extrapolate 
the high density evolution, but not as much as the Zel'dovich 
approximation, where caustic formation tends towards a too 
high distribution. In contrast to both models the Edgeworth 
approximation underpredicts the high density regions and 
develops a large wiggle in the positive density regime. 



Figure 4. Comparison with other distributions. We normalise all 
the models to the lognormal distribution (finorm(<5); dotted line), 
since we find that this best fits our N-body simulation results over 
a wide range of parameter values. The plot shows the ratio (P(<5) — 
Plnorm('5))/^'inorm ('5) for our second-order calculation (solid line), 
the Zel'dovich approximation (Kofman et al. 1994; dot-dashed 
line) and the Edgeworth expansion (dotted line). The variance 
used for the comparison was (tq = 0.25. 



1995), and that the lognormal makes a very useful fitting 
function. 

Our second-order PDF fits the peak fairly accurately, 
but tends to be slightly too skewed. This effect becomes 
greater at higher variance. The tails of the distribution 
match the N-body simulations rather well, although at high 
variances and low densities the distribution is too high. This 
is due to the real distribution being constrained to go to 
zero at zero density, whereas, as we have already remarked 
the second-order PDF is not. In comparison the Edgeworth 
expansion, taken to second order (Juszkiewicz et al. 1995, 
Bernardeau & Kofman 1995) 



1 



27r 



(31) 



where Hnix) is a Hermite polynomial, is also slightly too 
skewed and tends to undershoot the high density tail com- 
pared to the simulations. At large variance the expansion 
breaks down and the Edgeworth distribution develops a wig- 
gle on the high density tail. In fact the Edgeworth expansion 
fares badly even when taken to third order because the un- 
physical wiggles are exaggerated and the low density tail 
becomes more rapidly negative. 



4 COMPARISON WITH OTHER 
APPROXIMATIONS 

Finally we compare our PDF with other approximations 
for the 1-point density distribution. In figure ^ we show 
three distributions, our second-order PDF calculation (solid 
line), the Zel'dovich approximation (Kofman et al. 1994; 
dot-dashed line) and the Edgeworth expansion (dotted line) 
plotted relative to the lognormal distribution. Our choice of 
the lognormal distribution as the point of normalisation is 



5 SUMMARY 

We have presented a new method for calculating the non- 
linear evolution of the 1-point density distribution function 
using the Chapman-Kolmogorov equation to propagate the 
probability distribution, and second-order perturbation the- 
ory to evolve the density field. This has the advantage over 
other methods that the resultant probability distribution 
must be positive definite, and can be readily extended to in- 
clude the effects of Eulerian deterministic or stochastic bias 
and redshift space distortions. The main disadvantages of 
our method are that it is not obvious how to include the 
effects of smoothing in the final distribution and it is diffi- 
cult to see how one could extend the evolution of the den- 
sity field to higher order. However, despite these problems 
we find extremely good agreement with numerical simula- 
tions and that our distribution is an improvement on other 
approximations. Since we derive the characteristic function 
the moments and cumulants (connected moments) of the 
field are easily derived, as is the distribution of any lo- 
cal transformations of the density field. The formalism we 
have introduced can also be used to calculate more compli- 
cated 1- and 2-point distributions of cosmologically interest- 
ing fields. 
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